library(sf) # Manejo y análisis de datos espaciales en formato vectorial
library(terra) # Procesamiento y análisis de datos raster
library(gstat) # Modelado geoestadístico e interpolación espacial
library(tidyverse) # Conjunto de paquetes para manipulación, visualización y análisis de datos
library(tmap) # Producción cartogràfica en formato visor y plot
library(spatialEco) # Analisis espacial avanzadoProbabilidad de quema
103039 - Gestió de Riscos en la Planificació Forestal - UdL
1 Antecedentes
El incendio forestal del Pont de Vilomara, ocurrido en julio de 2022, fue uno de los más devastadores de la temporada en Catalunya, arrasando más de 1.700 hectáreas de terreno en la comarca del Bages. Este incendio se propagó con gran rapidez debido a las altas temperaturas, la sequía extrema y la sequedad del combustible, poniendo en peligro numerosas viviendas y obligando a la evacuación de cientos de personas.
Una de las zonas más afectadas fue la urbanización River Park, donde el fuego avanzó sin control y consumió decenas de casas, dejando a muchas familias sin hogar. La intensa labor de los bomberos permitió contener las llamas antes de que la destrucción fuera aún mayor, pero los daños materiales y el impacto fue significativo. River Park se convirtió en un símbolo de la vulnerabilidad de las urbanizaciones ubicadas en la Interfase Urbano Forestal ante el riesgo creciente de incendios forestales, generando un debate sobre la necesidad de mejorar las medidas de prevención y gestión del territorio.
Por ello, este trabajo se centra en la simulación de la longitud de llama, por ende la energía liberada, y en la propuesta de tratamientos que permitan reducirla en las inmediaciones de las zonas habitadas. A través de este ejercicio, buscamos identificar medidas efectivas de gestión del combustible y planificación del territorio que minimicen la propagación del fuego y protejan a las comunidades ante futuros incendios.
2 Archivo Landscape
El proyecto PREVINCAT es un proyecto liderado por el Centro de Ciencia y Tecnología Forestal de Cataluña (CTFC), junto con los Bomberos de la Generalitat, para proporcionar información básica para la modelización de incendios forestales y planificar futuras actuaciones para la prevención de incendios forestales.
Un archivo de paisaje (.LCP) consiste en un ráster en formato multi-banda, que se utiliza para las simulaciones del comportamiento y efectos de los incendios forestales. Concretamente la información que debe contener es:
Elevación
Pendiente
Orientación
Modelos de combustible (Scott and Burgan, 2005 - Tabla en Campus Virtual)
Cobertura arbórea
Altura del arbolado
Altura de la primera rama
Densidad aparente de las copas
Variables que componen el fichero Landscape (.lcp) Fuente: Finney, 2005
2.1 Definición del área de interes (AOI)
El proyecto PREVINCAT nos permite acceder a dicha información para toda Catalunya. PREVINCAT utiliza las Zonas Homogéneas de Régimen de Incendio (ZHR) para distribuir la información. Estas ZHR son descargables en el siguiente enlace y debéis descargar las que intersecten con nuestra área de estudio: El área quemada en el IF Pont de Vilomara + un buffer de 2km.
Para ello, primero de todo cargaremos las librerías necesarias para trabajar. Recordad que si es la primera vez que utilizáis la librería debéis instalarla primero utilizando la función install.packages().
2.2 Acceso y mosaicado información previncat
En segundo lugar cargaremos nuestra área de estudio, el Incendio Forestal de el Pont de Vilomara, disponible en Campus Virutal, le aplicaremos un buffer de 2km y lo definiremos como AOI. Se aplica dicha distancia para evitar el efecto borde que puede producirse en los limites del AOI.
AOI <- st_read("./Imputs/Perimetre_PontVilomara_25831.shp",quiet=T) %>% st_buffer(2000)En tercera instancia, se debe cargar y analizar que ZHR intersectan con nuestro AOI. Para ello, primeramente se carga el shapefile de con las delimitaciones espaciales de la ZHR utilizando la función st_read() y posteriormente se evalúa la intersección con la función st_filter().
ZHR <- st_read("./Preparacion/ZHR2014/ZHR.shp", quiet=T)
ZHR_sel <- ZHR %>% st_filter(AOI)
message(paste("ZHR seleccionadas:",paste(ZHR_sel$ZHR, collapse = ",")))ZHR seleccionadas: 26,36,37tm_shape(ZHR)+
tm_polygons()+
tm_shape(ZHR_sel)+
tm_polygons(col="red")+
tm_shape(AOI %>% st_cast("LINESTRING")) +
tm_lines(col="blue")
Como se puede observar el área de estudio intersecciona con las ZHR 26,36,37 y debemos decargar todos los rasters de las 6 vairables mencionadas anteriormente. Una vez descargados se procederá al mosaicado y recorte por el AOI.
Advertencia
Los archivos deben guardarse y descomprimirse en una única carpeta
#Ruta donde se almacenan los archicos descargados y descomprimidos
path <- "./Previncat/"
# Creación de vector con los nombres de las variables
vars <- list.files(path,
pattern = ".asc$",
full.names = F) %>% substr(.,1,str_locate(.,"_")[,1]-1) %>%
unique(.)
# Creación de vector con las rutas a los archivos
files <- list.files(path,
pattern = ".asc$",
full.names = T)
#Mosaicar y recortar
for(i in vars){
r.files <- files[files %>% str_detect(i)]
r <- lapply(r.files,rast)
m <- mosaic(sprc(r)) %>% crop(.,AOI)
plot(m, main=i)
names(m) <- i
m <- project(m, "epsg:25831","near")
writeRaster(m,paste0("./LCP_vars/",i,".tif"),overwrite=T)
}
3 Puntos de ignición simulados
En el presente ejercicio simularemos 10.000 incendios iniciados en los puntos con mayor probabilidad de inicarse un incendio según el histórico de incendios.
Para ello debemos:
Seleccionar los puntos de ignición ocurridos en un radio de 10km del incendio
Estimar el umbral de distancia a partir del cual los incendios son diferentes en términos de área quemada
Calcular la densidad kernel normalizada de igniciones con una distancia de búsqueda establecida en el punto 2. Además se debe establecer como densidad 0 aquellos espacios no combustibles.
Generar una muestra aleatoria balanceada espacialmente según la densidad observada en el punto 3.
Info
Encontraréis los puntos de ignición registrados en la EGIF (1998-2015) en el apartado Recursos de Campus virtual
3.1 Puntos de ignición ocurridos en un radio de 10km
Para ampliar el rango espacial de análisis y tener más muestra que las posibles igniciones ocurridas en el AOI, se amplía la búsqueda de igniciones ocurridas dentro de un radio de 10km.
buff <- AOI %>% st_buffer(8000) # Previamente habiamos aplicado 2km al general el AOI
ignit <- st_read("./Imputs/Ignitions/Ignicions.shp", quiet=T)
ignit_20 <- ignit %>% st_filter(buff)
tm_shape(ignit)+
tm_dots("grey")+
tm_shape(ignit_20)+
tm_dots("brown")
#Eliminar de la memoria todas las igniciones
rm(ignit)
# Vaciado de la memoria de papelera
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 3579957 191.2 6996292 373.7 6996292 373.7
Vcells 5680972 43.4 12345926 94.2 12345925 94.23.2 Dependencia espacial del área quemada
El área quemada es un indicador de la forma en la que queman los incendios y de las características topográficas, climatológicas y de vegetación que rigen en una zona. Es por ello que se busca ese umbral de distancia a partir del cual los incendios muestran un comportamiento diferente. Ese valor será usado como rádio de búsqueda en la función de densidad desarrollada en el siguiente punto.
Para ello, primero se calibra el semivariograma. Donde se representa la semivarianza a lo largo del gradiente de distancia. Esta, por lo general, tiende a crecer hasta que se estabiliza. El punto de estabilización nos marca ese umbral de distancia a partir del cual los incendios son diferentes en términos de área quemada y todo las variables que subyacen bajo ella. Se calcula el semivariograma empírico que es derivado de muestras observaciones, de la siguiente manera:
Se establece una relación o formula
variable~1.La fuente de información
cutoffhace referencia a la distancia máxima a evaluarwidtha los cortes cada \(x\) unicades de distancia
ve <- variogram(SizeHa~1,ignit_20,cutoff=10000, width=250)
head(ve) np dist gamma dir.hor dir.ver id
1 107 133.9027 19.560791 0 0 var1
2 142 375.0521 7.853788 0 0 var1
3 200 625.5520 395.953401 0 0 var1
4 217 877.3185 697.608888 0 0 var1
5 314 1129.0535 1508.768882 0 0 var1
6 342 1375.3032 2207.342806 0 0 var1plot(ve, plot.numbers = T)
Una vez tenemos el semivariograma empírico, ajustamos el semivariograma teórico. Que representa una curva con una ecuación conocida. Debemos seleccionar aquella curva que se ajuste mejor a la distribución de puntos observada en el semivariograma empírico:
En este caso en particular, el semivariograma spherical (Sph) es el que ajusta mejor
vt <- fit.variogram(ve, vgm("Sph"))
plot(ve, pl = T, model = vt)
kableExtra::kable(vt%>% select(model,range))| model | range |
|---|---|
| Nug | 0.000 |
| Sph | 7494.745 |
Como se observa en la tabla anterior, la distancia a la que los incendios empiezan a presentar diferencias en aria quemada es a los ~7,5km.
3.3 Densidad de puntos
Finalmente se deben generar una densidad de puntos donde se debe establecer un rádio de búsqueda máximo de 7,5km.
# Establecer como referencia espacial un raster cualquiera del archivo LCP
r.ref <- rast("./Preparacion/LCP_vars/cbd.tif")
# Cálculo de la densidad kernel
dens <- sf.kde(
x = ignit_20,
ref = r.ref,
res = 20,
bw = vt$range[2],
standardize = T,
scale.factor = 1
) %>% resample(r.ref,"near")
UBmask <- rast("./Preparacion/LCP_vars/model.tif")
UBmask[UBmask<=100] <- 0
UBmask[UBmask!=0] <- 1
plot(UBmask)
dens <- dens*UBmask
names(dens) <- "dens"
varnames(dens) <- "kernel density"
crs(dens) <- crs(r.ref)
plot(dens)
plot(AOI %>% st_geometry(),add=T, border="red")
3.4 Generación de muestra aleatoria espacialmente balanceada
Creación de una muestra aleatoria balanceada por la densidad. Las especificaciones de FlamMap requieren de los puntos de ignición simulados tengan las siguientes características:
nombre de las columnas
id, X_coord, Y_coordseparador
,decimal
.
sample <- dens %>% as.data.frame(xy=T) %>%
slice_sample(n=10000, weight_by = dens) %>%
select(-dens) %>%
mutate(ID=1:nrow(.) %>% as.integer()) %>%
rename(X_coord=x,
Y_coord=y) %>%
relocate(ID,X_coord,Y_coord)write.table(sample,"./Preparacion/pt10000.txt",sep = ",", dec=".",row.names = F)
Adjuntando el paquete: 'tidyterra'The following object is masked from 'package:stats':
filter
4 Simulación del riesgo de incendio en FlamMap
4.1 Generación del archivo landscape (*.lcp)
Ejecutar en FlamMap
Landscape > Create Landscape file…
Importante
Revisar unidades en las que PREVINCAT sirve los datos y establecerlas correctamente en el desplegable
4.2 Definición de características ambientales
Los vientos y la humedad del combustible son los principales factores que propician la propagación de incendios forestales. Dado que queremos evaluar una situación extrema, fijaremos el percentil 97 de las condiciones de viento y humedad observadas como umbral. Este valor representa un escenario de alto riesgo, pero no es un caso completamente atípico. En este contexto, se elige porque permite simular situaciones de peligro elevado pero realista, proporcionando una base sólida para la planificación y gestión de incendios forestales.
Según los datos aportados por la Agencia Meteorológica Meteocat, la velocidad de viento dominante es N
Por otro lado, debemos extrapolar el valor del P97 según la tabla de % de datos de vv que aparece en el margen inferior de la imagen de la rosa de los vientos mediante el siguiente procedimiento
tableVV <- data.frame(inf=c(0,.5,2.5,5,7.5,10,20),
sup=c(0.5,2.5,5,7.5,10,20,Inf),
freq_datos=c(18.8,62.4,16.3,2.5,0,0,0)) %>%
mutate(frec_accum=cumsum(freq_datos))
print(tableVV) inf sup freq_datos frec_accum
1 0.0 0.5 18.8 18.8
2 0.5 2.5 62.4 81.2
3 2.5 5.0 16.3 97.5
4 5.0 7.5 2.5 100.0
5 7.5 10.0 0.0 100.0
6 10.0 20.0 0.0 100.0
7 20.0 Inf 0.0 100.0\(P_k = L_i + \left( \frac{k - F_{i-1}}{F_i - F_{i-1}} \right) \times (L_s - L_i)\)
\(P_k\) valor de la variable en el percentil deseado.
\(L_i\) es el límite inferior del intervalo donde se encuentra \(P_k\).
\(L_s\) es el límite superior del intervalo.
\(F_{i-1}\) es la frecuencia acumulada antes del intervalo.
\(F_i\) es la frecuencia acumulada al final del intervalo.
\(k\) es el percentil deseado (en este caso, 97).
calcular_percentil <- function(tabla, percentil) {
# Encontrar el intervalo donde se encuentra el percentil
for (i in 1:nrow(tabla)) {
if (tabla$frec_accum[i] >= percentil) {
# Extraer valores del intervalo correspondiente
Li <- tabla$inf[i] # Límite inferior
Ls <- tabla$sup[i] # Límite superior
Fi_1 <- ifelse(i == 1, 0, tabla$frec_accum[i - 1]) # Frecuencia acumulada anterior
Fi <- tabla$frec_accum[i] # Frecuencia acumulada actual
# Aplicar la fórmula de interpolación lineal
Pk <- Li + ((percentil - Fi_1) / (Fi - Fi_1)) * (Ls - Li)
return(Pk)
}
}
return(NA) # Devuelve NA si el percentil no se encuentra en la tabla
}
VVP97 <- calcular_percentil(tableVV,97)
print(VVP97)[1] 4.923313Por tanto la velocidad extrema de viento se sitúa a los \(4.92 m/s \space (17.712 km/h)\) y la dirección de viento dominante es Norte.
4.2.1 Determinar niveles de humedad del combustible
Por otro lado, debe establecerse el contenido de humedad en la vegetación. Este se calcula mediante indices y ecuaciones que contemplan trayectorias temporales de humedad relativa, temperatura y precipitación.
En este caso utilizaremos los percentiles de humedad del combustible estimados para la Vegueria de la Catalunya Central por Alcasena, 2019.
| FuelModel | x1h | x10h | x100h | LW | LH | VV | DV |
|---|---|---|---|---|---|---|---|
| 0 | 8 | 9 | 13 | 50 | 60 | 4.92 | N |
5 Determinar la duración de los incendios simulados
Para intentar recrear al máximo la realidad de la zona de estudio, se buscará cuan es el tiempo optimo en alcanzar el tamaño medio de los incendios de la zona.
Por ello se buscan los incendios que hemos determinado anteriormente que eran parecidos en area quemada, aquellos que se encuentran a ~7,5 Km del AOI y calculamos la media de sus superficies quemadas.
ignit_20 %>% st_filter(AOI %>% st_buffer(vt$range[2])) %>%
st_drop_geometry() %>%
filter(SizeHa>10) %>%
summarise(mean(SizeHa)) %>% pull()[1] 101.3135Por tanto debemos recrear 10.000 incendios de unas que de promedio quemen 101 ha con las condiciones ambientales extremas establecidas establecidas. Para ello iremos incrementando el tiempo de quema y hasta alcanzar el umbral de tiempo necesario para quemar esa misma área promedio.
Ejecutar en FlamMap
Analysis Area > New FlamMap/MTT/TOM run con los parámetros
Debéis ajustar todos los parámetros requeridos (marcados en rojo)
Finalmente, debéis repitiendo simulaciones e ir incrementando el tiempo hasta alcanzar una área media quemada similar al histórico.
Hay ejecutaremos el algoritmo Minimum Travel Time (MTT), que, de manera simplificada, calcula el comportamiento del fuego buscando el conjunto de caminos con tiempos mínimos de propagación del fuego a partir de igniciones proporcionadas (ampliar información). Estableciendo el tiempo de quema calibrado y marcando las casillas para seleccionar los outputs deseados
Outputs a seleccionar
Burn probabilities
Perímeters
FLP-20 Bin Metric
Fire Size List
Muy importante
Se deben exportar los archivos des de el panel lateral de FlamMap y abrirlos con QGIS o ArcGIS, para los objetos espaciales y Excel o R para las tablas. No se recomienda abrirlos en el propio software ya que al ser pesado el programa puede no responder y perderse todo el trabajo realizado.
Parte 2: En esta segunda parte se procesan, tratan e interpretan los resultados obtenidos.
6 Outputs MTT
En el siguiente punto, se listan los productos Los productos que nos devuelve el algoritmo MTT:
Perimetros MTT: Los 10000 perímetros de las diferentes igniciones simuladasBurn Probabilities: Ráster con la probabilidad
FLP metric: Longitud de llama condicional. Probabilidad de tener una longitud de llama determinada.
Fire Size List: Tabla con las coordenadas de cada ignición y el tamaño de incendio en cada una de ellas.
| FIRE_NUM | XStart | YStart | FireSize.Acres | FireSize.Ha |
|---|---|---|---|---|
| 0 | 410595.5 | 4621331 | 362.9484 | 146.88013 |
| 1 | 409436.7 | 4613521 | 236.7270 | 95.80008 |
| 2 | 410649.4 | 4621957 | 537.8002 | 217.64019 |
| 3 | 410258.6 | 4624639 | 594.9309 | 240.76021 |
| 4 | 406216.4 | 4614296 | 196.2017 | 79.40007 |
| 5 | 404922.9 | 4616681 | 323.5104 | 130.92012 |
La longitud de llama responde a parámetros físicos como la intensidad de la llama (kW/m) que está ligada al calor de combustión o a la tasa de consumo de la vegetación. Se estima a partir de la expresión propuesta por Byram (1959):
\[ L_f=0.0775·I_b^{0.46} \]
Donde:
\(L_f\) es la longitud de llama
\(I_b\) es la intensidad de Byram
6.1 Cálculo longitud de llama
# Leer los datos de Probabilidad de Longitud de Llama (FLP) y ordenarlos en orden descendente por PBurn
LFP <- read_csv("./Outputs/FLP.csv") %>% arrange(desc(PBurn))Rows: 256055 Columns: 23
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (23): XPos, YPos, PBurn, FIL1, FIL2, FIL3, FIL4, FIL5, FIL6, FIL7, FIL8,...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# Crear una secuencia de constantes para el cálculo de la longitud de llama
constantes <- c(seq(0.5,9.5,.5)-0.25, 12)
# Calcular la longitud de llama (FL) y convertir los datos a un objeto espacial 'sf'
FL <- LFP %>% select(XPos,YPos,PBurn) %>% bind_cols(
sweep(LFP[, 4:ncol(LFP)], 2, constantes, `*`) %>% mutate(FL=rowSums(.)) %>% select(FL)) %>%
st_as_sf(coords = c("XPos","YPos")) %>%
st_set_crs(25831) # Establecer el sistema de coordenadas
# Leer el raster de referencia y ajustar su sistema de coordenadas
ras.ref <- r.ref
crs(ras.ref) <- "epsg:25831"
# Rasterizar los datos de longitud de llama (FL) utilizando el raster de referencia
FL <- rasterize(FL, ras.ref, field="FL") %>% mask(.,AOI)
# Visualizar el raster resultante
plot(FL %>% crop(AOI))
class_mat <- matrix(c(-Inf,1,1,
1,2.5,2,
2.5,3.5,3,
3.5,Inf,4),
ncol=3, byrow = T)
FL_clas <- classify(FL,class_mat)
coltab(FL_clas) <- data.frame(values = c(1,2,3,4),
cols = c("darkgreen", "gold", "orange", "brown"))
mylevels <-
levels(FL_clas) <- data.frame(values = c(1,2,3,4),
cats = c("Bajo", "Medio", "Alto", "Muy Alto"))
plot(FL_clas)
RiverPark<- st_read("./Imputs/RiverPark.shp", quiet=T)
RiverPark_fill <- st_multipolygon(lapply(st_geometry(RiverPark) , function(x) x[1]))
plot(FL_clas %>% crop(RiverPark %>% st_buffer(200)))
plot(RiverPark_fill %>% st_geometry(), , add=T, border="blue", lwd=2)
6.2 Tratamientos
Para reducir la exposición de bienes a altas intensidades de fuego se suelen realizar una serie de tratamientos forestales, como por ejemplo:
Reducción de combustibles de superficie
Quemas prescritas. Bosques de coníferas y pastos arbustivos. Requiere tecnificación, riesgo de escape. Ventana meteorológica. Plan de quemas.
Tratamientos mecánicos. Interfaz urbano-forestal, bosques de frondosas y zonas de elevada carga de combustibles. No requiere tecnificación.
Tratamiento de fracción arbolada
Aplicación de poda baja.
Eliminación de la tangencia en copas.
Conservación de áreas tratadas
Uso de ganadería extensiva en sierras o montes cercados. Grandes superficies.
Uso de pastoreo dirigido para la conservación de cortafuegos o puntos estratégicos
En este caso concreto estos tratamientos se podrían plasmar como:
Alterar el Modelo de combustible, por ejemplo de SB4 a SB1
Reducir la FCC hasta el 50%
Subir la base de la copa hasta \(3 m\)
Reducir la densidad aparente de copas hasta \(0.08 kg/m^3\)
Reducir la altura del arbolado mediante claras selectivas.
7 Ejercicio a realizar
Evaluar cual es la longitud de llama en las proximidades de la urbanización (un buffer de 60 m 3 píxeles).
Explorar cuales son los valores de las variables utilizadas en el archivo
*.lcpen las inmediaciones de la urbanización (buffer 60m).Diseñar tratamientos para reducir la longitud de llama en las inmediaciones de la urbanización.
Replantear los rasters del archivo landscape y reflejando los tratamientos.
Correr de nuevo la simulación de 10000 igniciones con FamMap incorporando los nuevos rásters